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Spatially-distributed, nonequilibrium chemical systems described by a Markov chain model are 
considered. The evolution of such systems arises from a combination of local birth-death reactive 
events and random walks executed by the particles on a lattice. The parameter 7, the ratio of 
characteristic time scales of reaction and diffusion, is used to gauge the relative contributions of 
these two processes to the overall dynamics. For the case of relatively fast diffusion, i.e. 7 < 1, 
an approximate solution to the Markov chain in the form of a perturbation expansion in powers 
of 7 is derived. Kinetic equations for the average concentrations that follow from the solution 
differ from the mass-action law and contain memory terms. For a reaction-diffusion system with 
Willamowski-Rossler reaction mechanism, we further derive the following two results: a) in the 
limit of 7 — » 0, these memory terms vanish and the mass-action law is recovered; b) the memory 
kernel is found to assume a simple exponential form. A comparison with numerical results from 
lattice gas automaton simulations is also carried out. 
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I. INTRODUCTION 



The dynamics of chemical reactions in condensed me- 
dia is commonly described by a set of reaction-diffusion 
equationsa of the following general form: 



<9x(r, t) 
dt 



R(x(r, f)) + DV 2 x(r,f), 



(1) 



where x(r, t) is vector of concentrations, r position vec- 
tor, t the time; R(x(r, t)) denotes mass-action law terms, 
and D is the matrix of diffusion coefficients. 

Such a description, while able to capture some of the 
features of the system's evolution, is not entirely free 
of shortcomings. For instance, as can be easily seen 
from ([!]), rescaling the diffusion coefficients by a con- 
stant factor is equivalent to rescaling of length. Contrary 
to this observation, the reaction dynamics in spatially- 
distributed systems varies with diffusion in a much more 
complex way; thus, it has been shown that in the limit of 
slow diffusionj-jjeaction dynamics has pronounced mem- 
ory character pa while in the case of very fast diffusion it 
closely follows the mass-action law. Significant progress 
has been made in understanding the dynamics of simple 
diffusion-limited chemical reactions (e.g. A + B — > 0) 
in low-dimensional systems.u Also, Chapman-Enskog 
methods have been applied to multi-variate master equa- 
tions to study reactive correlations in simple reaction- 
diffusion systems.EJ However, to date the dynamics of 
more complex reaction-diffusion systems whose mass- 
action laws exhibit periodic and chaotic oscillations has 
received much less attention. Recent numerical studies 
on such systems usiug lattice gas automaton modelso and 
master equationaJ'El have indicated that there are inter- 
esting effects pertaining to modification of the structure 
of oscillatory and chaotic attractors with decreasing rate 
of diffusion. 



The purpose of this paper is to develop a statistical 
mechanical theory that approximately accounts for the 
influence of finite diffusion and reactive correlations on 
the dynamics of nonequilibrium oscillatory chemical sys- 
tems. We use a discrete-space, discrete-time probabilistic 
model defined on a set of sites, or nodes of a lattice. Rs 
evolution is due to a combination of local reactive events 
and random walks which are executed by the particles of 
reactants on the lattice. 

Our approach is based on a Chapman-Enskog-like de- 
velopments applied to a Markov chain model. The ex- 
pansion parameter 7 is the ratio of characteristic time 
scales of diffusion and reaction processes, small for the 
diffusion-dominated regimes. We show that the corre- 
sponding kinetic equations for spatially-averaged concen- 
trations involve memory terms. As an example, we ap- 
ply thiSj-theory to Willamowski-Rossler reaction-diffusion 
system.E3 We find that the memory kernel for this 
system assumes simple exponential form. The mem- 
ory kernel should have the same functional form for any 
reaction-diffusion system whose reaction mechanism in- 
volves mono- and bimolecular steps only. The dynam- 
ics predicted by the kinetic equations indicates that the 
influence of diffusion essentially amounts to a backward 
shift in the bifurcation cascade of the mass-action law. 
We compare these theoretical predictions with lattice gas 
automaton simulations of Willamowski-Rossler model, 
and find that they are qualitatively in agreement. 

The paper is organized as follows. In Sec. |l| we intro- 
duce the Markov chain model and construct its evolution 
operator. In Sec. Ill we present a perturbation expan- 



sion for the probability distribution function around the 
local binomial form and obtain equations that determine 
evolution of the terms of this expansion. The perturbed 
distribution function incorporates the local correlations 
between different species induced by the reactive events 
and thus cannot be factorized into a product of single- 



1 



species distributions. The derivation of kinetic equations 
for global average concentrations (i.e. local concentra- 
tions averaged over the volume of the system) concludes 
this section. In Sec. IV we derive the memory kernel for 
Willamowski-Rossler reaction-diffusion system and com- 
pare the predictions of the theory with the results of lat- 
tice gas automaton simulations. Finally, Sec. ^ contains 
discussion of our results. 



II. MARKOV CHAIN DYNAMICS 

The system under consideration consists of a number 
of chemical species diffusing and reacting in solution. We 
assume a discrete-space, discrete-time description of the 
system by partitioning the space into cells of volume V 
and time into intervals of arbitrary length. In the follow- 
ing, we view each cell as a node on a lattice and operate 



within this lattice representation. The solvent comprises 
a chemically inert species so that its only effect is to 
maintain the random walk dynamics of the solute par- 
ticles. Furthermore, we associate N distinct "channels" 
with every cell; each channel can be occupied by no more 
then one particle of each solute species. Thus, if the num- 
ber of solute species is v, the number of solute particles 
in each cell cannot exceed vN - this is the so-called ex- 
clusion principle. At any given moment of time, the state 
of the system is fully specified by a set of state variables - 
the number and chemical nature of the particles occupy- 
ing every cell. The system evolves by transition to one of 
the accessible states once every time interval with proba- 
bility determined by the initial state. From the statistical 
point of view the system is described by the probability 
distribution function P(s, t) defined on the space of dis- 
crete states of the system, whose time dependence is due 
to a Markov chain Jld i.e.: 



1) - P(s, n) = YyT ss ,P{s', n) - T s , s P(s, »)], 



(2) 



where s is a set of state variables, n is the integer time 
variable, and the summation extends over all possible 
states. Clearly, due to discreteness of the state variables, 
eq. (0) can be also written in matrix form: 



P s (n + 1) - P s (n) = J2 W ss <P s >{r 



(3) 



where W ss i are the matrix elements of evolution opera- 
tor, defined in terms of transition probabilities per unit 
time as follows: 



T 

1 s« 



(4) 



Here S s r s is a Kronecker delta. 

Henceforth, all vector quantities, such as state vari- 
ables or concentrations, are denoted by bold lowercase 
letters, while bold uppercase letters stand for matrix 
quantities. Components of vectors and matrices are spec- 
ified where necessary by lower indices. We designate op- 
erators by a hat over corresponding symbols, and aver- 
age quantities by an overbar. Angle brackets will often 
be used as a shorthand notation for summing over the 
ensemble of state variables. 

For reaction-diffusion systems, transitions between dif- 
ferent states occur due to two competing processes: pas- 
sive transport of particles between cells and reactive col- 
lisions which change the nature of the particles within 
cells according to the reaction mechanism. Correspond- 
ingly, the evolution operator for the present model can 
be written as a sum of two terms, one for the diffusion 
process and the other for reactions. Below we derive the 
particular form of reaction and diffusion evolution oper- 
ators to be used in our study. 



A. Diffusion Markov chain 



Let the particles at every node be randomly distributed 
among ./V channels. We pose a diffusion rule accord- 
ing to which the particles that have identical chemical 
nature and belong to the same channel propagate syn- 
chronously and in the same direction between neighbor- 
ing nodes. The choice of a particular channel, direction 
of propagation, and species to be propagated is made 
at random each time the diffusion evolution operator is 
applied. We also impose periodic boundary conditions 
which, together with the above, completely characterizes 
transport of matter within the system. 

A similar diffusion model was recently used in sim- 
ulations of FitzHugh-Nagumo system.E3 It was shown 
that the occupancies of individual nodes become statis- 
tically independent of each other after a short period of 
relaxation. Hence, for long-time behaviour, correlations 
between nodes are negligible, and the full probability dis- 
tribution can be factorized into a product of single-site 
distributions. This reduction of description is readily ex- 
tended to the present model; specifically, one has for non- 
vanishing reduced transition probabilities: 



rpD = x i( r ) (, _ Xi(r,n) 

Xi(r)-l,tfi(r) vN I jy 



D 



vN 



I - 



N 



(5) 



where 
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r'eAf(r) x(r') 

= — V Si(r',n), 

r'eAf(r) 

to is the coordination number of the lattice and 7V(r) is 
neighbourhood of r. Thus, Xi( r T n ) is an average of the 
mean occupation number of the i-th species over the im- 
mediate neighbourhood of node r. The factors in the first 
expression in (|5|) simply reflect the fact that the proba- 
bility of a diffusive hop of a particle of the i-th species 
from the node r to a neighbouring node is given by a 
product of the probability that the i-th species is chosen 



for the hop, l/V, the probability that a chosen channel 
is occupied at the node r, Xi(r)/N, and the probability 
that the same channel is empty at the neighbouring node, 
1 — Xi( r i n ) /N. A similar interpretation can be also given 
to the second expression in (||). 

The reduced probability distribution P(x(r),i) evolves 
according to the following equation: 

P(x(r), n + 1) - P(x(r), n) = W D P(x(r), n) 

= E W D * ( II % r)) x,(r) I P W r )< n), (6) 

« \j¥* / 

with matrix elements of W D given by (cf. (H) and (0)): 



W. 



D' 

Xi (r) ,x'. (r) 



vN 



N 
N 



vN 



vN 
N 



1 



N 



^x'. (r),Xi (r) + l 



(7) 



x 1 - (r),Xi (r) • 



Thus, the structure of diffusion Markov chain (||) arises 
naturally from the diffusion rule posed in the beginning, 
the existence of a maximum occupation number at any 
node (i.e. the exclusion principle), and the reduction of 
description to the level of single-node distribution func- 
tions. 



Kinetic equations for average occupancies (concentra- 
tions) can be obtained in the standard way, by multi- 
plying both sides of (^|) by Xi(r) and summing over all 
possible states. One finds, after a few simple transforma- 
tions: 



Xi(r, n + 1) - Xj(r,7i) 



1 



[xi(r',n) - Xi(r,n)} 



'QV(r) 

D Axi(r, n), 



(8) 



where A is the discrete Laplacian operator. This is just 
a discretized diffusion equation with diffusion coefficient 
D = (vmN)^ 1 (the same for all species). 



details of the derivation £ar the models with exclusion 
principle can be found in;H3 here we only briefly summa- 
rize the general ideas and give the final results. 



B. Reaction Markov chain 

The probabilities of reactive transitions are defined 
similarly to those in lattice gas automaton models. The 



We consider reactions which are strictly local, i.e. they 
involve only particles occupying the same node. Let the 
overall reaction mechanism consist of r elementary steps 
of the form: 



+ n{ 2 X 2 + ■■■ + ntjXj 



nl-iXi + nUXa 



where the indices i — 1 , . . . , r and j — 1 , . . • , v la- 
bel steps and species, respectively. The transitions due 
to the i-th step in the mechanism are characterized by 



two sets of stoichiometric coefficients {r 



,,n; 



} and 



,,n 



'„}. Specifically, if we define the vector m 



(i) 



with elements 



the forward reaction at 



3 



node r can be represented in terms of occupancies by: 

x(r) — > x(r) + mW. 
Similarly, for the reverse reaction: 



x(r) — > x(r) — m 



(i) 



The transition probabilities of these birth-death reactive 
events can be derived using a simple combinatorial ar- 
gument. That is, one expects that the probability of 
reaction must be proportional to the number of ways in 
which the particles of the reactants can be combined in 
order for that reaction to occur, given the present state 
of the node. Specifically, we put: 



J x(r)+m( i ),x(r) 



T 



R 



x(r)-mW ; x ( r ) 



V 



N\ 

{N-vQ\ Xj (r)\ 



(xj(r) - u 



- n f . )! 

i,3 J 



(9) 



The prefactors in (||) are chosen so that the mean- 
field kinetic equations can be written in a neat ana- 
lytic form. The exclusion principle forbids those reactive 
events which produce particles of any species in excess of 
N. Therefore, we set: 



^x(r)+m(*> ,x(r) — ^) if Xj 



mf > N, 



T. 



0. if Xj 



3 



> N 



(10) 



x(r)-m< i ) ,x(r) " "3 ""3 

for any i,j. Substituting the transition probabilities (|^), 
modified according to , into (Q) yields the matrix el- 
ements of reaction evolution operator W . 



Let us consider the mean-field dynamics of the re- 
action Markov chain. To begin, we note that in the 
mean-field limit (i.e. in the limit of a well-stirred sys- 
tem where the rate of diffusion is infinite), the particles 
are re-distributed among the nodes of the lattice instan- 
taneously at every moment of time. Hence, the local 
correlations between different species are effectively de- 
stroyed at the same moment as the local reactive events 
create them. This implies that the particles populate 
every node independently of each other, and the local 
probability distribution is the time-dependent multivari- 
ate binomial: 



P(x(r),n) = P B (x(r),n) = J] ' ^ 



N 



N 



N-Xj(r) 



(ID 



We have dropped the notation due to spatial coordinates that the kinetic equations in this limit have the following 
in the average concentrations in (|ll]) since the system is form: 
homogeneous in the mean-field picture. It is easy to show 



j(n+ l)-%(n) =^2xjW R P B (x,n) 



i=l 

^■(x(n)). 



h3 



h n x k i,k ( n ) - k -* n x k i ' k ( n ) 



fe=i 



k=l 



+ Q i (x(n)) 



(12) 



One observes that the first term on the second line of (|12| ) 
is the mass-action law (in the discrete-time notation) ap- 
propriate for the overall reaction mechanism. The ad- 
ditional term Q(x(n)) describes the deviations of the 
mean-field dynamics from the mass-action law due to the 



exclusion corrections to the birth-death transition prob- 
abilities (^) . Eq. (|l^) is the exact mean- field rate law for 
our model. In the application to Willamowski-Rossler 
reaction diffusion system in Sec. IV we shall consider 
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conditions where Q(x(n) is negligible and eq. ( jl^ ) re- 
duces to the mass- action law. 

For a finite diffusion rate, neither flll] ) nor (||) will 
be valid. The nature of this breakdown of mean-field 
description forms the central point of the present study. 
Note also that the binomial distribution ([ll]) depends on 
time only through its means Xj (n) which evolve according 
to (|l2[). We will use this time-dependent binomial dis- 
tribution as the zeroth-order approximation to the full 
distribution in the following section. 



III. PERTURBATION THEORY 



the value of diffusion coefficient D. Such a regime can be 
characterized by two time scales, tjj and tr, associated 
with diffusion and reaction processes, respectively. We 
consider finite systems so that long wavelength diffusion 
modes are suppressed and we may insure that td *C tr. 
In the context of our model td can be the mean time 
required for a particle to travel the distance between two 
neighbouring nodes, and tr the inverse of the smallest 
eigenvalue of reaction evolution operator W R . The devi- 
ation from the mean-field limit is measured by the ratio 
^P-, which will be denoted as 7 henceforth; 7 is equal to 
zero in the mean-field limit, and increases away from it. 



Suppose that the system is not completely homoge- 
neous, i.e. diffusion occurs at a finite rate determined by 



If 7 is small, the local probability distribution 
P(x(r),i) obeys 



P(x(r), n + 1) - P(x(r), n) = (~/W R + W D ^ P(x(r), n), 



(13) 



where W R and W D are the reaction and diffusion evo- 
lution operators constructed in previous section. We ex- 
pect that the solution of (|l^) monotonously approaches 
the local binomial distribution ([ll]) of the mean-field 
limit as 7 tends to zero. For a non-zero 7 this mean- 
field result is not valid, and corrections to the local bi- 



nomial distribution must be introduced. In view of this, 
for small 7 reactions can be considered as a small per- 
turbation to the pure diffusion process and perturbation 
theory can be applied to solve equation ([i"3|). Thus, we 
write the local probability distribution in the form of a 
perturbation series: 



P(x(r), n) = P s (x(r), n) + 7 i\(x(r),n) + 7^ 2 (x(r), n) 



(14) 



with the local binomial distribution as the zeroth-order gives: 
term. Substituting (Jl4|) into evolution equation ( |i~3| ) 



Pn(x(r),n + 1) - P s (x(r),n) + ]T 7 fc {P fc (x(r), n + 1) - P fe (x(r),n)} 

fc=i 

/ 00 
= (7^ + ^) P B (x(r),n)+^7 fe P fc (x(r)/, 



(15) 



fc=i 



in (|j~5| ) are equivalent to a discrete derivative of the bino- 
In order to extract the time dependence of each of mial distribution Ps(x(r),n) with respect to time. Re- 
the expansion terms from (|l5|), we use a perturbation calling that this distribution depends on time implicitly, 
theoretic, .procedure similar to Chapman-Enskog devel- through its means x(r, n), one can re- write the derivative 
opment.Bli-1 Specifically, we note that the first two terms as follows: 



Pg(x(r), n + 1) — Ps(x(r), n) = x(r, n + 1) — x(r, n 



— P B (x(r),n). 



(16) 
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The quantity x(r, n + 1) — x(r, n) appearing in 
is given by the kinetic equations corresponding to (|15[). 
The latter can be obtained in the standard way, i.e. mul- 



where K(x(r, n)) stands for mean-field dynamics (cf. eq. 
©), and 

(x(r)W^P fe (x(r), n)) - ^ x(r)W R P fe (x(r), n). (18) 

x(r) 

Before we proceed with the development of the solu- 
tion to the evolution equation ( |l5| ) one important ques- 
tion must be addressed, specifically that of the validity 
of the regular perturbation expansion (|lj) for the case of 
far-from-equilibrium chemical systems. It is well-known 
that the dynamics of such systems may be of relaxational 
character, i.e. it exhibits well-defined fast and slow parts. 
This is manifested by the significant jumps in reaction 
rates as well as the higher derivatives of concentrations 
with respect to time which are observed as the system 
switches from fast to slow dynamics, or vice versa. In 
principle, it is quite possible that such jumps can cause a 
similar singular behaviour in the probability distribution 



YPi(x(r),n+ 1) - 7 
+ 7 



and for Pi(x(r),n), 1 < i < k (cancelling out the factor 



The solution of this system of inhomogeneous linear equa- 
tions will be unique if we require: 

(p(x(r),n)) =0, (22) 
<x(r)Pi(x(r),n)) =0, i=l,...,k. 

Together with the solvability conditions (p2|), ( po|) and 
( pl[ ) form a closed hierarchy from which one may evalu- 
ate the terms of the perturbation expansion ( |l4| ) to any 
desired order. 



tiplying both sides of (|15| ) by x(r) and summing over all 
possible combinations of the occupation numbers. Using 
the properties of operators W D and W R introduced in 
Sec. 0, we find: 



(17) 



P(x(r),n), thereby making the description in terms of 
a regular perturbation series invalid. In order to avoid 
such breakdown, we require that 7 be small enough for 
the following relations to hold at all r and n: 

P(x(r),n) « 7 -\ (19) 
(x(r)IU K P(x(r),n)) i = 1,2,3,... 

With this requirement satisfied, the regular perturbation 
method remains valid, at least on the level of description 
provided by equations (|l^) and (p7|). 

Now, suppose that we want to determine the probabil- 
ity distribution up to the fc-th order term. For that pur- 
pose, it is sufficient to retain in ( p"7| ) terms up to the order 
7 fe . Then, with the use of Jl7j ) and (|l6|), we can separate 
( |l5| ) into a set of k equations, each describing evolution of 
one of the expansion terms Pi(x(r),n), . . . ,Pfc(x(r),n). 
Thus, the equation for Pi(x(r),n) reads: 



P B (x(r),n) (20) 
7 W D P 1 (x(r),n), 



7'): 



(21) 



In the remainder of this paper our attention will be 
limited to the first two terms of the series (|l4|). There- 
fore, we need only solve the lowest-order equation of 
the above hierarchy, i.e. (|20|). By direct computa- 
tion of M /£> PB(x(r) i n) we may show that the first term 
on the r.h.s. of ( |20| ) is zero and, consequently, that 
Pi(x(r),n + 1) — Pi(x(r),n) does not contain contribu- 
tions of order unity. For simplicity, we assume that the 
system consists of a single species. Using the definition 



00 

x(r, n + 1) - x(r, n) = D Ax(r, n) + 7K(x(r, n)) + 7 fc+1 ^x(r)TU i? P fe (x(r), n) 

k=l 



Pi(x(r),n) = 



W — D Ax(r, n) — 
ox 



W R - K(x(r,n)) A 
ox 



P B (x(r),n) + 



P(x(r),n + 1) - P(x(r),n) = TU fl p_i(x(r), n) + W D P t (x(r),n) 

- (x(r)iy fl P(x(r),n)) Ap s ( x (r),n). 



G 



of the diffusion evolution operator W D (see @ and (0)) we can write this term as: 



X(r,n) 



N 



1 - 



x(r) — 1 \ , , . . ( x(r) , 

1 - V L I Pb{x{v) - l,n) - 1 - 4/ P B (x(r),n) 



X(r,») 
iV 



TV 
a;(r) + 1 



TV 



(23) 



a;(r) 



Y PB(x(r)+l,n)-^P B (x(r),n) 



%(r,n) - x(r,n) 9 
TV 9x 



P B (x(r,n). 



Taking advantage of two easily proven results, 

1 — I Ps(a;-l,n) = V AT _, \ JJ P B (x,n), 



N 



Nx(n) 



(24) 



d xP B (x,n) - (x+l)P B (x+l,n) 

— P B {x,n) = — — . 

ox x(n) 



we can show that 



vanishes. Generalization of this 



result for a many-species system is straightforward due 
to the properties of the diffusion evolution operator (cf. 
(§))• 

With the first term on the r.h.s. of ( pfj| ) eliminated, we 
can obtain Pi(x(r),n) by solving a simple initial- value 
problem. Suppose that at n = the local probability dis- 
tribution is purely binomial, i.e. Pi(x(r),0) = 0. Then, 
one finds from (20): 



n-l 

P 1 (x(r),n)= (l + W D ) 

n'=0 



w 



K(x(r,n')) 



<9x 



Pfl(x(r),n'), 



(25) 



where I is the identity operator (I x > x — <5 X 'x), and 

^ \ n—n — 1 

I + W D ) is to be understood as a time-ordered 



product of n — n' — 1 terms ( I + W D ^j , with the leftmost 

term taken at the moment of time n — 1 and the right- 
most at n' + 1. In Appendix A we prove that (^H|) does 
satisfy the solvability conditions (p2]). 

The quantities whose dynamics will be of interest to 
us are the global concentrations x(t), i.e. the local con- 
centrations x(r, t) averaged over all nodes of the lattice. 
(Henceforth, whenever we omit r in the arguments of 



the concentrations, the global concentrations are meant.) 
The kinetic equations for the local concentrations are 
given by ( |l7j ) ; at the level of present approximation we 
can truncate the infinite series on the r.h.s. of ( |l7| ) after 
the term of order 7 2 . To obtain kinetic equations for the 
global concentrations we substitute Pi(x(r), n) from (Efj) 
into (fl7"l), sum both sides over all r, and divide by the to- 
tal number of nodes M. Since diffusion conserves the 
total number of particles in the system, the pure diffu- 
sion term in (|l^) yields zero on summation and we obtain 
for the global concentration of the i-th species: 



i(n + 1) - 3?i(n) = ^Y, K ^ ™)) + ^ E E (xi(*)W R (i 



W 



n — n' — 1 



SP B ( X (r),ri] 



(26) 



where S = W R — K(x(r, n'))-J=, and x(n) is the vector put: 
of global concentrations. 

If diffusion is sufficiently fast (7 <C 1), the system re- 
mains nearly homogeneous at all times, i.e. Xi(r,n) are 
slowly varying functions of r which do not deviate sig- 
nificantly from their spatial averages Xi(n). Hence, if we 



Xi(r,n) « Xi(n) (27) 

for all r and n, the summation over r in ( p6| ) becomes 
trivial, and we find: 



n-l 

Xi(n+1)-Xi{n) = 1 K l {5c{n)) + 1 2 ^ ( Xi (r)W a (l + W 

n'=0 



n—n' — 1 



(28) 
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where the local concentrations appearing implicitly in 
W D , S, and Pb( x :^) m the last term are replaced by 
global concentrations x(n). 

The kinetic equations (Pq ) do not readily offer a deeper 
insight into dynamics of the model since the memory 
function in the last term cannot be easily evaluated for 
general reaction mechanisms. However, as we will show 
in the following section, for specific reaction mechanisms 
it is possible to write these equations in a tidy analytic 
form which allows clear physical interpretation. 

IV. APPLICATION TO 
WILLAMOWSKI-ROSSLER MODEL 

We now apply the above formalism to a reaction- 
diffusion system with Willamowski-Rossler reaction 



mechanism. This mechanism consists of the following 
elementary steps: 

Ai + X ^ 2X, X + Y ^ 2Y, 

fe_i k-2 

A 5 +Y ^ A 2 , X + Z 4 A 3 , 

K — 3 k — 4 

Ai + Z ^ 2Z . 

fe-5 

Here A\, A±, A§ are initiators, and A 2 , A3 are final prod- 
ucts; concentrations of all these species are held constant 
by external fluxes. The intermediates whose dynamics is 
followed are X, Y, and Z. The mass-action law for this 
model has the following form (indices 1,2,3 refer to X, 
Y , and Z species, respectively): 



— = KlXl{t) - K-lx\(t) - K 2 X 1 (t)x 2 (t) + K- 2 x\{t) - K4,Xl(t)x S (t) + K_4 = #1 (#1 (t) , X 2 (t) , X 3 (t)) , 

—— = K 2 xi(t)x 2 (t) - n- 2 xj(t) - K 3 X 2 {t) + k_ 3 = R 2 {xi(t),x 2 {t)) , (29) 
at 

— = -K 4 Xi(t)x 3 (t) + K 5 X 3 (t) - K- 5 xl(t) + K_4 = R3(Xi(t),X 3 (t)) , 



where the concentrations of initiators and products 
are incorporated in the values of rate constants Kj, 
i = 1, . . . , 5. Equations ( p9|) have been considered in 
earlier studies and are known to give rise tOiSHseriod- 
doubling cascade leading to chaotic attractorBt£ltL3 

We begin our analysis by calculating the contributions 
that the reactive steps of Willamowski-Rossler mecha- 
nism make to the memory function in (p8|). Below, the 
contributions from mono- and bimolecular steps are con- 
sidered in turn. 



A. Monomolecular steps 

All monomolecular steps in Willamowski-Rossler 
model have the following form: 

k 

Ai —> rijXj + mXi + n m A m , n,j,ni,n m = 0, 1. 

The matrix elements of the reaction evolution opera- 
tor corresponding to such steps can be obtained from 
the general expression for reactive transition probabili- 
ties (Q). We find: 

-hK (l - rrij6 Xj<N ) (1 - miS xltN ) , if x' = x, 
W^t = { Jik, if x' = x — m, 
0, for all other x', 

where the Kronecker deltas account for the exclusion 
principle, the concentration of species A4 is incorporated 



in the value of k and vector m is defined through the 
stoichiometric coefficients as outlined in Sec. ||. 

At this point, in order to facilitate further calcula- 
tions, we neglect the terms in the memory function that 
arise from the exclusion principle. Such approximation 
is justified for Willamowski-Rossler model, since those 
terms can always be made arbitrarily small by a simple 
rescaling of concentrations. As an example, suppose that 
iV = 10 and the concentrations are scaled so that they 
practically never exceed 1.5 particles per node. Given 
that the reactive transition probabilities are of the or- 
der of 10 -3 - 10 -4 , numerical estimates show that the 
terms due to the exclusion principle are 2-3 orders of 
magnitude smaller than all the other terms in (|28| ) . (The 
concentration scaling and the value of N used in the lat- 
tice gas automaton simulations of Willamowski-Rossler 
model below are the same as considered here.) 

The contribution to the memory function can now be 
evaluated in a fairly straightforward way. We substitute 
the evolution operator from above in place of the leftmost 
W R in the last term of (|28|) , and calculate the average. 
After a few transformations, we arrive at: 



7 hK 



n-l 

E 

n'=0 



I + W 



D 



(30) 



One immediately recognizes in (|3^) one of the two quanti- 
ties (to within a constant multiplier) that we have dealt 
with in a different context in Appendix A (cf. (A.l)). 
Using the results obtained in Appendix A we find that 
( p^p| ) vanishes. Therefore, we conclude that monomolecu- 
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lar reactive steps in Willamowski-Rossler mechanism do 
not contribute to the memory function in 



—hkxiXj (l — m ' 2 +1 S Xij N) > if x' = x, 
Wxx' — { hk (xi — mi) (xj — rrij) , if x' = x — m, 
0, for all other x', 



B. Bimolecular steps 

The bimolecular steps in Willamowski-Rossler model 
fall into three classes; these are listed below, along with 
the corresponding matrix elements of the W R operator: 

k 

1) 2Xi—*Xi+rijXj+TnAi, nj,m — 0, 1, 



3) Ai + Xj A2Xj, 



W R , 



-hnxj (l — S Xjt jA , if x' = x, 



Jin (xj 
0, 



if x' = x — m, 
for all other x', 



y y xx' 



N 

y-i 

N S\hk (xi — m,i) (x,i — mi — 1) , if x' = x — m, 



jt-lhkXi (xi - 1) (l - m,jS x ^ N ) , if x' = x, 
JV 



0, for all other x', 



2) Xi + Xj A mXi + njAj, m = 0, 2; nj = 0, 1, 



where rrij = 1 and the concentration of Ai is incorporated 
into the value of k. 

We neglect the terms arising from the exclusion prin- 
ciple on the basis of the rescaling argument given above, 
and calculate the contributions to the memory function 
due to these reactive steps. We find: 



N 



n — n' — l 



SP B (x,ri) 



2) j 2 hk (xiXj (i + w D 

n'=0 ^ 
n— 1 / 

3) j 2 hn (xj^i + V 



n — n — 1 



SP B (x,ri) 



(31) 



71 — 71 — 1 



SP B (x,n')). 



Note that the last contribution in (|3l]) is identical (to 
within a constant multiplier) t o on e of the quantities 
considered in Appendix A (cf. (A..9)). From the results 
of Appendix A we infer that this contribution vanishes. 



Also, using the method of Appendix A one can show that 
the first two contributions in ( |3l| ) can be reduced to the 
following memory terms with a simple exponential ker- 
nel: 



i\t ™ — 1 / i \ n—n'—l 



2 > ^El 1 -^ 



i>N 



n — n' — l 



x 2 SP B {x, n'h)), 



XiXjSPB{~x-, n'h) 



(32) 



Summarizing the results obtained above, we conclude 
that the only reactive steps that contribute to the mem- 



ory function in (28) are those involving reactions between 
the particles of intermediate species. The contribution is 
given by the first term in ( |32| ) if the reacting particles 
are of the same chemical nature, and by the second term 

otherwise two ^ or ^ ne ^ species: 



2X ^ 


A x - 




2Y ^ 


X + 


Y, 


X + Y H 


2Y, 




x + z % 


A a , 





By inspection of the mechanism, we find four con- 
tributing steps for the X species: 



2Y -4 x + y, 
X + Y h 2Y, 
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and two for the Z species: 



2Z k ^ A 4 + Z, 



X + Z 



-4, 



To obtain the kinetic equations for the global concentra- 
tions of each species, one simply adds contributions ap- 
propriate for the steps in each of the three groups above 
and substitutes the result in place of the memory term 
in (28). Proceeding in this way we arrive at: 



xi(n + 1) - xi(n) = jhRi(x(n)) - ^hjJ-^-j ^ (1 - A)" ™ 1 n-iC Xl , Xl (n') - n~ 2 C X2 . X2 (n') 



n'=0 



n-l 



x 2 (n + 1) - x 2 (n) = jhR 2 {x(n)) - j 2 h ]T (1 - A) n_w _1 

n-l 

x 3 (n + 1) - x 3 (n) = 7 Wk(x(n)) - £ (1 - A)™ - ™' -1 



^ 2 h^2 (1 - \) n n 1 K 2 C XuX2 {n') + KiC Xl , X3 (n') 

Cx 2 ,x 2 ( n ) ~~ K 2 C XliX2 (n ) 



n'=0 
n-l 



N - 1 



NK-5 

N-l 



C X3 , X3 (n') + K i C Xl ^ x . i {n') 



(33) 



(34) 
(35) 



where 



S=W R - fcR(5E(n)) 



d_ 



(36) 



and Rj (x(n)), i = 1,2,3 are the mass- action law terms 
from (E9L) , A = = ^ for Willamowski-Rossler model. 
The procedure which we used above to evaluate the mem- 
ory terms in (Pq) by neglecting terms due to the exclusion 
principle can be extended to any reaction mechanism in- 
volving only mono- and bimolecular steps. Note that in 



the scaling limit where the terms due to the exclusion 
principle can be omitted in evaluating all average quan- 
tities, the mean-field reaction dynamics is entirely due to 
the mass-action law, i.e. K(x(n)) = hR(x(n)) (cf. (p"2|)). 
The continuous-time form of eqns ( |33"| ) - ([55]) is derived 
in Appendix B. 

In order to evaluate the averages appearing in the 
equations ( J33| ) - ( |35| ) one needs to know the matrix el- 
ements of the full reaction evolution operator W R . The 
latter are defined through the reactive transition prob- 
abilities which, for Willamowski-Rossler model, are as 
follows: 



1 xi + l,X2,x 3 ;x 



hnxxi (1 - 6 Xu n) ■ 



T 



N 



T 



X\— 1,X2 + 1,X3 ;x 



flK 2 XiX2 (1 - 5 X2jN ) , T R + 1 



N-l 
N 



x 2 -l,x 3 :x ]\f _ I 



hn^iXi (xi — 1) , 

hn- 2 x 2 (x 2 - 1) (1 - S XltN ) . 



rpR 

X\ ,X2 — I.X3 ;x 

rpR 



T 



x\ — l,x 2 ,X3 — l;x 
R _ 

Xl,X2,X3 + l,X 



= h,K 3 x 2 , 

/IK4X1X3, 



Tx-l^ + I.xs-^ - hK -3 (1 - &x 2 ,n) 



rpR. 

± x± + l,x 2 , a; 3 + l;x 



hn 5 x 3 (1 - 6 X3 , N ) . 



T 



hn- A (1 - 5 Xl . N ) (1 - 5 X3 ,n) ■ 
N 



x,i,x% ,X3 — l;x 



N - 1 



/1K-5X3 (x 3 - 1) 



Here x = (xi, x 2 , x 3 ). Using the transition probabilities 
([37]) to compute C XuXj , we obtain: 



(37) 
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C XuXl (t) = 2h 
C X2l x 2 \t) = 2ft 
Cx 3 ,x 3 (t) = 2ft 

C" Xl ,x 3 (t) — h 



iV 



Ea(*) 



K 2 a;i(t)x2(t) - ( 1 - -rr ) + -rr^W 



a; 3 (t)\ k_ 4 . 



iV 



K_ 2 si(t) (1 



TV / AT 
2z 2 _ , . m A x 2 (t) -Zi(i) 



(38) 



iV 



k_ 4 - K 4 .xi(t)a;3(t) 1 



- K 2 xi(t)x 2 (t) ( 1 
A" 



N 



Equations (|33|) - (|35[) together with fl38|) describe the dy- 
namics of the global concentrations within terms of order 
7 2 for Willamowski-Rossler reaction-diffusion system. 

We conclude this discussion by proving that the dy- 
namics predicted by eqns ( |33"| ) - ( ]3q ) is consistent with 
the initial assumptions of the theory, i.e. that it con- 
verges to mass-action law dynamics in the well-stirred 
regime (i.e. 7 -> 0). Because ft is arbitrary, we observe 



that assigning a smaller value to 7 simply amounts to 
following the dynamics due to ([33]) - (^5|) on a smaller 
time scale. This implies that under a rescaling of time 
ft' = 7ft, the limit 7 — > is formally equivalent to the 
limit ft' — > with 7 fixed. Defining a new discrete (real- 
valued) time variable t„ = -ft', we can re-cast ( J33] ) - ( |35| ) 
in the following form (we put 7 = 1 and consider only 
eqn. (|34|), for simplicity): 



n-l 



J'2 



(t„ + ft') -S 2 (f») = ft'-R 2 (x(t„)) - ft' 2 J] (1 - A)" 



n'=0 



iV «_ 2 C X2>X3 (f 



N - 1 



Keeping i ra and t n i constant as we let ft' tend to zero, we obtain: 



(39) 



dx 2 (t) 
dt 



i? 2 (x(i)) - lim 

ft'— *0 



t — h 



N K_ 2 C X2 ^ X2 (t') K 2C xl , X2 (t') 



N - 1 



dt' 



(40) 



The result on the last line follows since the kernel of the 
integral term vanishes exponentially for all t — t'. Similar 
results hold for eqns |3^) and ( ^B|) as well. 



0.01, K5 = 16.5, k_5 = 0.5; this choice corresponds to 
the period-4 regime in the mass-action dynamics. The 
integration time step ft was equal to 5 • 10 . 



C. Numerical results 

We have solved numerically the continuous-time form 
of the equations (|3^) - ( |35| ) (cf. Appendix B) using Eu- 
ler's algorithm; the integral terms were evaluated at every 
time step by Simpson's method. In Fig. |l] we present the 
phase space trajectories obtained for different values of 
the diffusion coefficient (as measured by the exponent of 
the memory kernel A) as well as for the mean-field regime 
(A = 00). The values of rate constants used in these 
calculations are as follows: k\ — 31.2, K-\ — 0.2, k 2 — 
1.533, k_ 2 = 0.1, k 3 = 10.8, k„ 3 = 0.12, k 4 = 1-02, k_ 4 = 



It is evident from Fig. [l] that, as the diffusion becomes 
slower (i.e., the value of A becomes smaller), the trajec- 
tories are monotonously transformed so that they tra- 
verse the period-doubling cascade of the mass-action law 
( p9| ) in the backward direction: the smaller the value of 
the diffusion coefficient (with rate constants unchanged) , 
the less the dynamics of the global concentrations corre- 
sponds to the mass-action law. 
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A = 00 A = 0.5 




A = 2.0 A = 0.25 




FIG. 1. Phase space trajectories obtained from equations 
( p3[ ) - (|3^) for different values of A. Mass-action law trajec- 
tory (A = oo ) is period-4. Indices 1, 2, and 3 in the coordinate 
labels correspond to X, Y, and Z species, respectively. All 
trajectories are drawn to the same scale. The values of A are 
given in units of i (h = 5 • 1CF 4 ). 



where m^j, mu are positive integers, and / is the iden- 
tity operator. Eq. ( |4l| ) is approximately equivalent to the 
reaction-diffusion evolution equation (p"3|), with parame- 
ter 7 given, to the leading order, by the ratio 22s., This 
fact, along with definition of the parameter A (A = h 
for Willamowski-Rossler model), allows us to approxi- 
mately relate the results of an automaton simulation with 
a certain value of 11111 to a solution of the kinetic equa- 
tions (|3^) - (f35|) with an appropriate value of A. 

In Fig. ^ we present phase-space trajectories from typ- 
ical lattice gas automaton simulations for the parametric 
regimes where the mass-action law dynamics is period- 1 
and period-2. All rate constants except K2 are the same 
as in Fig. [IJ n 2 is equal to 1.4 and 1.5, respectively. For 
each regime, A was equal to 0.667, 13.333, and oo, the 
latter corresponding to the limit of well-stirred system. 
The simulations of well-stirred system were performed by 
re-seeding the nodes of the lattice according to binomial 
distribution at every time step. 

We observe that the trajectories from these simulations 
exhibit the same general trend as those computed using 
kinetic equations ( |33| ) - ( |35| ) in Fig. EI. 



Furthermore, this breakdown of mass-action descrip- 
tion can be modeled by a parametric shift in the mass- 
action law (along with rescaling of time and, possibly, 
concentrations) such that the dynamics proceeds from 
the period-4 regime (corresponding to the mean-field 
regime in Fig. yj) to period-2, period-1, and finally to 
stable focus. 



We now turn to the comparison of the numerical solu- 
tions of equations ( |33| ) - (|3f]) and the dynamics observed 
in the lattice gas automaton simulations of Willamowski- 
Rossler model. The lattice gas automaton used in these 
simulations was implemented in essentially the same way 
as in the earlier studiesp'E3 the only difference is in the 
nature of the diffusion rule, which we modify so that it 
corresponds to the diffusion evolution operator ((?]). The 
simulations were performed for a triangular lattice of size 
200 x 200 nodes, with time step h = 5 • 10~ 4 . The ex- 
clusion parameter N was equal to 10 (i.e. no more than 
10 particles of each species were allowed to occupy any 
node at any time). All concentrations were scaled by a 
factor of 40 so that they do not exceed 1.5 particles per 
node throughout the simulations. 

In a lattice gas automaton simulation, the local prob- 
ability distribution evolves in time according to the fol- 
lowing equation: 



(41) 




FIG. 2. Phase space trajectories from lattice gas automa- 
ton simulations for K2 = 1.5 (left) and K2 = 1.4 (right). The 
mass-action law dynamics is period-2 and period-1, respec- 
tively. For each parametric regime, the value of A (in units 
of i) are: 0.667 (top panel), 13.333 (middle panel), and oo, 
(bottom panel). All trajectories are drawn to the same scale. 
Concentration scaling parameter is 40. 

Namely, note that for the regime where the mass- 



P(x(r),n + 1) = ll+ W R j (l + W D ) P(x(r),n), 
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action law dynamics is period-2, the trajectory is gradu- 
ally transformed from noisy period-2 (A = oo) to noisy 
period-1 (A = 13.333) and finally to the stable focus 
(A = 0.667). In order to investigate this effect in more 
detail, we construct a Poincare section for the trajecto- 
ries from automaton simulation. The half-plane of sec- 
tion (denoted by V) is shown in Fig. 0, bottom left panel, 
and defined as follows: {P : x\ > 0.3; ~x~2 — 0.25; xj, > 0}. 

In Fig. |^ we show the normalized distribution of con- 
centration of X species on Poincare section V for dif- 
ferent values A. The value of K2 is 1.4; corresponding 
mass-action law dynamics is period-1. 




0.1 0.3 SP , 0.5 0.7 



FIG. 3. Distribution of concentration of X species on 
Poincare section V for different values of A. The mass-action 
law dynamics is period-1 («2 = 1.4). The graphs are num- 
bered according to the value of A (in units of i) as fol- 
lows: 1,A = 0.667; 2, A = 3.333; 3, A = 6.667; 4, A = 10.0; 
5, A = 13.333. Each distribution is computed from a set of 
2000 data points. 

One observes that the peak of the distribution shifts to- 
ward lower concentrations as A decreases, indicating that 
the trajectory shrinks as the rate diffusion is decreased. 
That this shrinking is due to a parametric shift and not to 
simple rescaling of volume can be seen from the fact that 
the shift of the peak vanishes as the diffusion becomes 
faster and the automaton trajectory approaches that of 
the mass-action law. Note that the width of the distribu- 
tion does not vanish as the rate of diffusion increases, but 
instead approaches a non-zero limiting value. This effect 
arises as a result of fluctuations which the global aver- 
age concentrations experience due to the finite number 
of particles involved in the simulations. 

Fig. ^presents similar distributions computed for K2 = 
1.5 (the mass-action law dynamics is period-2). The val- 
ues of A used in these calculations are the same as in 
Fig. g. 




0.1 0.3 0.5 0.7 0.9 1.1 

—sec 

X l 



FIG. 4. Distribution of concentration of X species on 
Poincare section V for different values of A. The mass-action 
law dynamics is period-2 (ac2 = 1.5). The graphs are num- 
bered according to the value of A as in Fig. H. Each distribu- 
tion is computed from a set of 2000 data points. 



Note that for the highest value of A (i.e. the fastest dif- 
fusion) the width of distribution in Fig. |^ is considerably 
greater than in Fig. |^. The additional spreading comes 
from the fact that the dynamics is already beyond the 
period-doubling bifucation at this value of A. However, 
the concentration fluctuations destroy the bimodality in 
distribution, and we are unable to distinguish the loops 
of period-2 trajectory even for relatively fast diffusion. 
Larger system sizes are necessary in order to resolve this 
structure. 

Finally, in Fig. [5] we show the average oscillation am- 
plitude for the concentration of Y species for a series 
of lattice gas automaton simulations. The values of the 
rate constants correspond to period-1 regime in the mass- 
action law dynamics and are the same for all simulations, 
while A is varied from 0.667 to 13.333. The results ob- 
tained from the kinetic equations ([53j) - ([55]) are shown 
by the solid line. The oscillation amplitude x^ mp is de- 
fined as the difference between the maximum and the 
minimum of Xi over one full cycle of trajectory around 
the unstable focus. Since the present theory is approx- 
imate, our primary aim here is to present a qualitative, 
rather than quantitative, comparison of theoretical data 
and simulations. Therefore, the exponent A is given in 
Fig. U in rescaled form, A resc = A ^ Aref . Here A rc f is the 
reference value of A, which is equal to 0.667 for the au- 
tomaton simulations, and 0.128 for the theoretical data; 
at these values of A the amplitude obtained from simu- 
lations and predicted theoretically match exactly. 
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FIG. 5. Oscillation amplitude for the concentration of Y 
species in period- 1 regime for different values of A obtained 
from the kinetic equations (^) - l^^j (solid line) and from 
lattice gas automaton simulations (open diamonds). 

One sees that, qualitatively, kinetic equations ( |33| ) 
( |35|) correctly describe the breakdown of the mass-action 
law description for the dynamics of global concentrations 
as the diffusion coefficient is decreased. 



V. CONCLUSIONS 

In this paper, the breakdown of the mean-field dy- 
namics in non-equilibrium reaction-diffusion systems is 
studied by the means of a Chapman-Enskog type ex- 
pansion using a probabilistic discrete-space, discrete-time 
model. It is shown that this breakdown is related to 
the memory effects in the evolution of average concen- 
trations that emerge for finite values of diffusion coeffi- 
cient. The memory kernel, although quite complicated in 
general case, was simplified for the case of Willamowski- 
Rossler reaction-diffusion system to a simple exponential 
e -A(t-t ^ wnere t ue exponent A is determined by the dif- 
fusion coefficient and t is the (continuous) time variable. 

A numerical analysis of theoretically derived kinetic 
equations with memory led to a further insight into 
the nature of the breakdown of mean-field dynamics. 
Namely, the effect of a decreasing diffusion coefficient 
on the reaction dynamics can be described by a back- 
ward shift in the bifurcation cascade of the mass-action 



law, i.e. the breakdown is essentially a parametric phe- 
nomenon. This conclusion is further supported by a com- 
parison to the results of lattice gas automaton simula- 
tions of Willamowski-Rossler model. 

The breakdown of mean-field description for the dy- 
namics of reaction-diffusion systems has been considered 
in a number of earlier studies. Thus, simulations of the 
reaction-diffusion master equation for Brusselator model 
showed that the deterministic limit cycle is gradually de- 
stroyed as the rate of diffusion is decreased and the inho- 
mogeneous dynamic modes become excited.Q Recently, 
Bussemaker and Brito used ring kinetic theory to ex- 
plain shrinking of the limit cycle fox Maginu model with 
attenuation in the rate of diffusion.EJ These authors were 
able to qualitatively reproduce the effect of breakdown of 
the mean-field kinetics observed in lattice gas automaton 
simulations by accounting for the reactive correlations in 
a phcnomcnological way. Because the probabilistic model 
used in these studies is similar to that we used here, it 
may be interesting to consider the relation among these 
approaches. 

There is a number of important questions which were 
not considered in this work. For instance, development of 
an analytic means for estimating the effective parametric 
shift exhibited by the solutions of kinetic equations j3^ ) 
- ( |35| ) as a function of diffusion coefficient would be use- 
ful. A different, no- less-challenging question concerns the 
extension of this theory to the level of local concentra- 
tions and the incorporation of long wavelength diffusion 
modes into the description of the dynamics. Such an ex- 
tension would clearly be expedient for studies of stability 
and evolution of spatio-temporal structures in reaction- 
diffusion systems. 
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APPENDIX A 



Substituting Pi(x(r),n) from (Ul 
first of the solvability conditions (p 



i into the r.h.s. of the 
3) gives: 



Pi(x(r),n)) = W(i + W 



n'=0 



W H ~K(x(r,n')) — 



ft(x(r),n'; 



(A.l) 
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To simplify this expression, note that by definition: 



W D \ de- 



P' =1 7r(p') 



(A.2) 



where p is any positive integer and JDwp') 

notes sum over all possible time-ordered products of p' 
operators W D selected from the full set of p such opera- 
tors. 



Given (A2), we can write (Al) as: 



P 1 (x(r),n))=^ 



n'=0 



W H -K{x{r, n')) 



+ E E(M 

'=1 7T(/) N 



<9x 



P fl (x(r),n'; 



W w -K(x(r, n')) 



9x 



ft(x(r),n'; 



(A.3) 



The first summand in ( |A.3| ) vanishes since for all n we 
have: 



W R P B (x(r),n)j =Q, 
P B (x(r),n)\ = 1. 



(A.4) 



The remainder of our analysis is easier to carry out in 
matrix form. We define: 



B 



(0 

x(r) 



(n)=(w D ) 1 SP B (x(v),n), 



where 



S = W R -K(x(r,n)) — , 



(A.5) 
(A.6) 



(w D ) { 



I. 



and re- write (A.3) in terms of vector elements 5^(n) 
as follows: 



n— 1 n— n' — 1 

P 1 (x(r),n))= £ E E (W'^J) (»'))• ( A - 7 ) 

n '=° i=1 7T(0 

Using the definition of the matrix elements of the diffu- 
sion evolution operator given by (^) and ([?]), we expand 
the leftmost W in (A. 7) and obtain: 



n—X n—n—1 



P 1 (x(r),n))=Y, E 

«-° '=! 7T(0 



E (w D B^<-> 



(A. 



Proceeding in the similar manner with the other condi- tion in (p2|), we find: 



n- 1 



n'=0 
n-1 



x(r)Pi(x(r), n ) = J] (x(r) 7 + W c SP s (x(r),n') 



n— n — 1 



(A.9) 



n—n' — l I . 

^|( X (r)ff B (x(r), n '))+ ^ Wx(r)(w D ) $P B (x(r),n') 

™'=0 I '=1 7t(Z) 



The summand on the second line of (A.9) again vanishes 
because of the following relations: 



Xl (r)W R P B {x(r),n)) = i^(x(r,n)), 

s i (r)P fl (x(r),n)) = Zi(r,n). (A.10) 
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Hence, one can again expand the leftmost operator W 
and re- write (A.9) as follows: 



n — l n—tl—X 



x(r)P 1 (x(r),n)) = E E E Uv)W D B^\n') 



n'=0 1=1 



(A.ll) 



n—l n — n — 1 



^E E T,{^n")-x(r)]B^(n')), 
n'=0 i=l 7r(/) 



where n" is the moment of time at which the leftmost 
operator W D was applied (0 < n" < n — n' — 1), Bv.(t) 



vector elements i ?^ r ^ (n) where possible, and repeat the 



is defined by (A. 5). 

One can now separate the leftmost operator W D from 



cal culat ions i n (A.ll ). Taking advantage of the results 
in ( A.4 ) and (|A.1C| ), we obtain: 



n — l n — n' — l 



-. ft I ft ft J. 

x(r)P 1 (x(r),n)) = -^ £ E (^(^") - z(r)] BfeV)) , 



AT 



=0 


=1 7T(Z) 


n-1 




n'=0 




n-1 


: ( , 

(n — n — 


n'=0 





n—7i' — 1 

E E 

n — 'n' — 1 



'= 2 7T(Z) 



(A.12) 



n-1 n-n' — 1 / , ,\ / , \ ( 

n'=0 1 = 1 V / \ / 

= 0. 



APPENDIX B 



It is convenient to consider the continuous-time limit 
of kinetic eqns ( j33|) - ( j35[ ) in rescaled time units, h' = jh. 
Thus, we define a discrete, real-valued time variable 
t n = -h', and re-cast (|||) - ( |35| ) in the following form: 



1G 



xi{t n H ) - 

7 



/i'i?i(x(i n )) 



iV 



AT 



n'=0 



(1-A) *7 



111 — '"/, fl ^ 7 K — iCxi ,x\ (in' ) K — 2Cx2,X2 (^rt') 



1 " J ^ ^7 



(B.l) 



2 2 (£n H ) - x 2 {t n ) 

7 



= h'R 2 (x(t n )) 



n-l 



h' 2 E (! - A ) 



t n — 4 t—h'/i 



n'=0 



N-l h'h 



h'h 



(B.2) 



#3(£n H ) - a; 3 (t„) 

7 



tiR 3 {x(t n )) 



h' 2 ^ (1 - A) 



£ n — t r—h/'y 



n'=0 



iV K_5Ca;3 ja ; 3 (i n ') K4C XljX3 (t n ') 



iV - 1 &77 



+ 



h'/7 



(B.3) 



We now have to evaluate the limit of (B.l) - (B.3) as 
h! tends to in such a way that the diffusion coefficient 
remains unchanged. Formally, we can do so by letting 
A = Ah! = 7A/1, and keeping A constant as h' — > 0. 
From the point of view of statistical theory, such an op- 



eration is equivalent to replacing the diffusion transition 
probabilities (|^) with appropriate transition rates (such 
that the resulting diffusion coefficient is identical in both 
cases) as the diffusion Markov chain (|^) is replaced with 
a continuous-time Markov process. Proceeding in this 
manner, we obtain: 



dxi(t) 
dt 



dx 2 (t) 
dt 

dt 



fli (*(*)) 



Mm) 

iZa(x(t)) 



N 



N - 1 



-A(t-t') 



K-lC Xl ,x 1 (t') — K-2Cx 2 ,x 2 (t') 



dt' 



(B.4) 



-A(t-t') 



,-A(t-t') 



-A(i-t') 



K 2C Xl ,x 2 

(t') + KiC XliXs (t') 

Nk 



dt' 



N - 1 



Nk-5_ 
N - 1 



Cx 2 ,x 2 (t) K 2^xi,x 2 {t ) 



^x^.x^it ) ~h ^4^a;i .X3 (£ ) 



(B.5) 
(B.6) 



,-(*) 



,(*) 



where Ca^fa) = 
defined by (pq) and (f 



s). Eqns (B.4) - (B.6) represent 
the continuous- time form of the kinetic equations ( |33| ) - 
(|35|). The parameter A is defined as A = = ^ h , 
where h is the time unit of diffusion Markov chain (@) . 
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